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Abstract 

A new representation of tlie exact time dependent solution of the discrete 
master equation is derived. This representation can be considered as contrac- 
tion of the path integral solution of Haken. It allows the calculation of the 
probability distribution of the occurence time for each path and is suitable as 
basis of new computational solution methods. 



Typeset using REVTgX 



1 



INTRODUCTION 



The master equation [|T]-^ plays a very important role in the description of stochasti- 
cally behaving systems. In statistical physics it arises in the MARKOFFian limit |Q of the 
generalized master equation 0. The generalized master equation is obtained as a special 
representation of the Nakajima-Zwanzig equation, which results by projection of the voN 
Neumann equation for the statistical operator on the relevant variables of the considered 
system . Applications of the master equation reach from nonequilibrium thermodynam- 
ics |]9HlI|, over chemistry [p|,p!2| and biology [l^ to the social sciences |1T4|-[I8[]. 



In most concrete master equation is only numerically soluble. For a small num- 

ber of possible system states computer programs for the solution of a linear system of first 
order ordinary differential equations can be used. However, in many situations Monte Carlo 
simulations have to be applied, yielding only approximate results. An alternative ap- 
proach is the path integral solution discussed below. It provides an exact expression for the 
time dependent solution even of the generalized master equation, for which very few solution 
methods exist up to now. This expression is a suitable starting point for the development of 
new computational solution methods. It also allows the calculation of the probability that a 
system takes a 'desired' or a 'catastrophic' path, which is very important for technical and 
other applications. 

THE PATH INTEGRAL SOLUTION 

Let p{x^ t) with < p{x, t) <\ and Y^xVi.^^ ^) = 1 denote the probability oi the considered 
system to be in state x at time t. Further, for x' ^ x let 

, /, ,x p{x' ,t + M\x,t) 
w[x \x; t) = jm — (1) 

be the transition rate (i.e., the transition probability p{x',t + At\x,t) per time unit At) 
from state x to state x'. Then, the master equation for the temporal development of the 
probability distribution of the system reads 



dp{x^t) ^ J2 [yj(^x\x';t)p{x',t) -w{x'\x;t)p{x,t)]. (2) 
Having the hmit At in mind, this equation can be written in the form 

p{x,t + At) = p{x,t) + At ^ w{x\x';t)p{x' ,t) — At ^ w{x'\x;t)p{x,t) 

= + At|x',t)p(x',t) (3) 

x' 

with 

p{x,t + At\x' ,t) = 6xx' + At ^ [w{x\x"]t)Sx'x"] — Sxx'At ^ w{x"\x;t) , (4) 

which is also called a propagator, and the Kronecker function 

1 if X = x' 

■= { (5) 
otherwise. 

The solution of the master equation (0) is, therefore. 



X]v_i a;jv-2 a;o U=l 



N 



J9(Xo,to), (6) 



where xat := x, := to + iAt, and At := (t — to)/N. This representation is a consequence 
of the Chapman-Kolmogorov equation and it is the basis of Haken's approximate 
path integral solution of the master equation |pO |. 



In the following, we will restrict our considerations to time independent transition rates 
w{x'\x; t) = w{x'\x) and derive a contracted form of expression (^. In order to achieve this 
simplification we utilize the fact that the system will normally not change the state Xi over 
many time periods At. That means, within the time interval t := t — t^ the system will 
usually have changed its state Xi a finite number of times only. Let xq,xi, . . . ,Xn be the 
sequence of states the system shows between times to and t. Then, 

Cn '■= Xq Xi ^ . . . ^ Xn {7) 

with Xn '■= X shall be called the path of the system, and 
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p{Cn, t) = p{xo ^ . . . ^ Xn,T) (8) 

shall denote the probability that the system takes the path C„ during the time interval r. 
This implies the relation 

oo 

p{x,t) =p{x,to + T) = J2J2p(^n,r) 

n=0 C„ 

oo 

= Y1 ■■■ Yl Pi^o ^ . . . ^ x„, r) . (9) 

The representation (P) of the solution of the master equation (|^) is the desired contraction 
of formula @ and has been suggested by Empacher [^. In order to derive the explicit 



form of the path probabilities p{Cn,T), we will extend an idea of Weidlich p2|. For this 
purpose we rewrite the master equation in vectorial form 

dp{t) 



dt 



{W - D)p{t) (10) 



with 



p{t):={p{l,t),...,p{x,t),...r (11) 



and the matrices 



W = (W,,,) := [w{x\x')] 
D = (d,,,) := {wj,,,) , (12) 



where 



Wxx = w{x\x) := 0, 

Wx := w{x'\x) . (13) 

X'iy^x) 



Wx is the total rate of transitions from state x to all other states. 
Applying the Laplace transformation 



f{u):= j dte-^'^'-'^^fit) (14) 
to 



with a function, vector or matrix / to equation (0) we obtain 

upiu) - pito) = (W - D)p{u) , (15) 

that means, 

p{u) = {uE-W + D)-'p{to), (16) 

where E = (^6xx'^ is the unity matrix and denotes the inverse of a matrix M. The 
desired solution of the master equation (p!0|) is, finally, found by inverse Laplace transfor- 
mation 

c+ioo 

f{u) = i^/j du e«(*-*")/(tz) (17) 

c— ioo 

of p{u). It yields 

—c+ioo 

pit) = — I du e-"(*-*«) (D-uE- Wy^pito) 
2ni J 

— c— ioo 
—c+ioo 

= — / due-''^'-'''^[E-{D-uE)-^W]-\D-uE)-^p{to). (18) 

— c— ioo 

The geometric series 

oo 

[E-{D- uE)-^W]'^ = J2[iD- uEy'wr (19) 

n=0 

converges, if the constant c is chosen sufficiently large. Taking into account 

[{D - uE)-'Wlx' = - ur%x.w{x"\x') = ^^^^ (20) 



Wt — u 



and (ITsI) we obtain 



E E{[P - uE)-'wr{D - uE)-'}x„x,p{xo, to) 

n=0 xo 

t E ••• E p(xo,to) (21) 

i=0 



with 



^xxo if n = 

w{Cn) = w{x, ^ . . . ^ := <! " ^ , , ^ , (22) 

_[J_ 'w[Xj\Xj-i) II n > 1. 

Now, evaluating (|18|) by means of (^) and comparing the result with (P) leads to the explicit 
expression for the path probabilities p(C„, r). It reads |[16[| 



— C+lOO 



p{Cn,T) = — du — p{xo,to) . (23) 

-c-ioo 



This important result is in accordance with a detailed heuristic derivation of Empacher 



By application of the residue theorem integral (^) can easily be evaluated. For example, if 
all total rates are identical (i.e., = w), we get 



whereas we obtain 



p{Cn, t) = —e-'^^wiCMxo, to) , (24) 



p{Cn, r) = ± f "'"^^(^n) ^^^^^ ^ (25^ 



i=0 



n (^^. -^^j 



i=o 



if the total rates Wr^ are all different from each other. Expressions (p^ and (^) are poly- 
nomially increasing functions in r for r ^ and exponentially decreasing functions for 
r ^ oo. 



Note that for some cases the 'path integral' (|^), (^) can be analytically evaluated pT 
Moreover, a similar relation exists for the generalized master equation 

dp{t) 



dt 



to 

with 



J dt' [W{t - t') - D{t - t')]p{t') (26) 



w(x\x'; t — t') li x' ^ X 
W,At-t')--={ (27) 
otherwise, 



since a Laplace transformation of (P6|) leads to |]5| 



up{u) - p{to) = [W{u) - D{u)]p{u) . (28) 

In ([2^ ) and (^) one has, therefore, only to replace w{x'\x) by w{x'\x; u) and Wx by Wx{u) := 
ti;(a;'|x; u). 

PATH OCCURENCE TIMES AND COMPUTATIONAL METHODS 



In numerical evaluations of formula (|^), one can, of course, not sum up an infinite number 
of different paths. One will rather restrict the summation to the paths which give the most 
important contributions. In order to find out which paths are negligible, we will calculate 
the probability distribution 

p(r|C„) :^ (29) 

of the occurence time r of a path C„ given that this path is traversed. Here, we have 
introduced the normalization factor 

oo 

p{Cn) := J dT p{Cn, t) . (30) 


We expect, that p(r|Cn) is maximal for a certain time f(C„). On the one hand, for r ^ t(C„) 
there is not enough time for the n transitions Xt^i — > Xj to occur. On the other hand, for 
T ^ f(C„) the likelihood of further transitions 

Xn a^n+l — ^ . . . — > Xn+l (31) 

with Z > 1 is great. In order to determine the mean (t)c„ of the path occurence time t and 
its variance 

Oc. := ((r - {r)cj)c. = {t\„ - {{r)cS (32) 

we need the quantities 

oo 

{t')c„ := J dr T'p{T\Cn) . (33) 



By means of the residue theorem one obtains K 



1 "V^ in 
1 [ . Vj[Lri 



{rlc = du 



w{Cn) k\ 

p{Cn) 27ri 



i=0 



p(xo,)f:o) 



°' ° p(C„) V^m/ {w^^ -u) 



From this we can derive 



u=0 



p{Cn) = piXo, to)w{Cn) H 



1=0 -"-t 



(34) 



(35) 



r 



and 



n 1 



(36) 



(37) 



In order to get an approximate solution of p{x, t) = p{x, to + t), we have only to sum up 
the path probabilities p(C„, r) of paths C^, for which 



\r- (r)cj < aJ9c„ 



(38) 



holds. The parameter a depends on the desired accuracy of the approximation. 

In the following the computational solution method shall be explained in more detail. 
Let us assume to have chosen a = 3 for the accuracy parameter. Then, about 99 percent 
of the distribution function p(r|C„) are found between ti := (r)c„ — ciyfdc^ and T2 : = 
{'^)c„ + This implies that for every given r the paths Cn which fulfil condition 

(p8|) will reconstruct about 99 percent of the distribution function p{x, to + 't) defined by 
(^). Therefore, one possible method of numerically calculating the distribution function 
p{x, to + r) is the following: 
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1. For every x set p{x,tQ + r) := and calculate Wx from w{x'\x) according to (|T3]). 

2. Let n = and generate all paths Cq := xq of length n = 0. Define the relevance status 
of each path Co by setting ■s(Co) := 1. Calculate {t)co '■= l/'U^xo ^Co '■= l/(u^xo)^ 
every path Cq. If condition ( P8| ) is fulfilled, then add exp(— Wa;or)p(xo, to) top(xo, ^0+^)- 

3. Let n := n + 1 and generate all paths C„ := C„_i — > a;„ for which the subpath C„ 
is relevant, that means, s{Cn) = 1. Calculate (t)c„ := (t)c„_i + and 6'c„ : = 
^c„_i + for every generated path C„. 

4. If condition (^) is satisfied, then calculate p{Cn, r) according to ( pSD (or the corre- 
sponding explicit formula of (p3D) and addp(C„, r) to ^o+t)- If r < (t)c„ — cy^^cT 
is fulfilled, then cancel the relevance of C„ by setting s(C„) := 0. 

5. If there are no relevant paths C„ with s(C„) = 1 left, then stop the procedure. Other- 
wise pro cede with 3. 

More refined, but more complicated algorithms can be developed for computers with little 
memory. 

CONCLUSIONS AND OUTLOOK 

We have derived a contracted and explicit path integral solution for the discrete master 
equation. The result is also apphcable to the generalized master equation, and it allows the 
development of numerical methods of solution. 

The advantage of the path integral solution with respect to other solution methods is 
that one can calculate the probability of a system to take 'desired' paths or 'catastrophic' 
paths consisting of sequences of selected states only. This is of interest for prognoses (e.g., 
szenario techniques) or the (technical) control of a system's temporal evolution, especially 
if the probability distribution develops several maxima. 
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